{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Compare FWI result with true model for the Overthrust model__\n",
    "\n",
    "Daniel Köhn \n",
    "Kiel, 16/07/2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Import Libraries__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "from matplotlib.colors import LightSource, Normalize\n",
    "from matplotlib.pyplot import gca\n",
    "from pylab import rcParams\n",
    "from matplotlib import rc\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib.ticker as mtick\n",
    "import scipy.ndimage.filters\n",
    "from scipy.ndimage.filters import gaussian_filter\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Activate different post-processing options__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "GAUSSIAN=1;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__FD grid dimensions__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DH = 0.2\n",
    "NX = 500\n",
    "NY = 304"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Wavefield clip value__ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "clip = 5e-2\n",
    "#clip = 5e1\n",
    "vpmin = 60.0\n",
    "vpmax = 800.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Define fonts__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "FSize = 25\n",
    "font = {'color':  'black',\n",
    "        'weight': 'normal',\n",
    "        'size': FSize}\n",
    "mpl.rc('xtick', labelsize=FSize) \n",
    "mpl.rc('ytick', labelsize=FSize) \n",
    "rcParams['figure.figsize'] = 16, 10.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Read FWI result and true model__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#f = open(\"../start/Kleinneudorf_fsurf_100_smooth.vs\")\n",
    "f = open(\"../FWI_results/25_11_2017_p100_elast/modelTest_vs_stage_4.bin\")\n",
    "data_type = np.dtype ('float32').newbyteorder ('<')\n",
    "mod_true = np.fromfile (f, dtype=data_type)\n",
    "mod_true = mod_true.reshape(NX,NY)\n",
    "mod_true = np.transpose(mod_true)\n",
    "mod_true = np.flipud(mod_true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = open(\"taper_p100.bin\")\n",
    "data_type = np.dtype ('float32').newbyteorder ('<')\n",
    "taper = np.fromfile (f, dtype=data_type)\n",
    "taper = taper.reshape(NX,NY)\n",
    "taper = np.transpose(taper)\n",
    "taper = np.flipud(taper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#f = open(\"23_11_2017_bp_40_50Hz_offset_20m_p100/jacobian_S_image.bin\")\n",
    "f = open(\"27_11_2017_bp_40_50Hz_offset_10m_p100_FWI/jacobian_S_image.bin\")\n",
    "data_type = np.dtype ('float32').newbyteorder ('<')\n",
    "RTM_TD = np.fromfile (f, dtype=data_type)\n",
    "RTM_TD = RTM_TD.reshape(NX,NY)\n",
    "RTM_TD = np.transpose(RTM_TD)\n",
    "RTM_TD = np.flipud(RTM_TD)\n",
    "RTM_TD = scipy.ndimage.filters.laplace(RTM_TD) # suppress low-wavenumber artifacts in image\n",
    "RTM_TD *= taper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Apply Gaussian filter__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if(GAUSSIAN==1):\n",
    "    RTM_TD = gaussian_filter(RTM_TD, sigma=[1,6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Define Axis__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.arange(0.0, DH*NX, DH)\n",
    "y = np.arange(0.0, DH*NY, DH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Scale RTM result with depth__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "RTM_scale = np.zeros((NX,NY))\n",
    "RTM_scale += np.flipud(y)**4\n",
    "RTM_TD*=RTM_scale.transpose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Define SubPlot__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def do_plot(n, model, cm, an, title, vpmin, vpmax):\n",
    "    \n",
    "    ax=plt.subplot(1, 1, n)\n",
    "    extent = [0.0,NX*DH,0.0,NY*DH]\n",
    "    #plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "    rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "    ## for Palatino and other serif fonts use:\n",
    "    #rc('font',**{'family':'serif','serif':['Palatino']})\n",
    "    #plt.rc('text', usetex=True)\n",
    "    rc('text', usetex=True)\n",
    "    \n",
    "    if(n==1):\n",
    "        im1 = plt.imshow(mod_true, cmap=plt.cm.jet, interpolation='nearest', extent=extent, aspect=1)\n",
    "        plt.hold(True)\n",
    "\n",
    "    im2 = plt.imshow(-RTM_TD, cmap=plt.cm.gray, alpha=.30, interpolation='bicubic',\n",
    "                 extent=extent, vmin=-clip, vmax=clip, aspect=1)\n",
    "    \n",
    "    a = gca()\n",
    "    a.set_xticklabels(a.get_xticks(), font)\n",
    "    a.set_yticklabels(a.get_yticks(), font)\n",
    "    a.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.0d'))\n",
    "    a.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.0d'))\n",
    "    #plt.axis('scaled')\n",
    "    plt.title(title, fontdict=font)\n",
    "    plt.ylabel('Depth [m]', fontdict=font)\n",
    "    plt.xlabel('Distance [m]', fontdict=font)\n",
    "    #ax.set_xticks([]) \n",
    "    plt.gca().invert_yaxis()\n",
    "    \n",
    "    # add annotation\n",
    "    #if n!=2:\n",
    "    #    plt.text(0.5, 4.2,an,fontdict=font,color='white',size=20)\n",
    "    \n",
    "    # fit and label colorbar\n",
    "    #divider = make_axes_locatable(ax)\n",
    "    #cax = divider.append_axes(\"right\", size=\"2.5%\", pad=0.05)\n",
    "    #cbar = plt.colorbar(im1, cax=cax)\n",
    "    #cbar.set_label(an, fontdict=font, labelpad=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Plot SubPlots__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.close('all')\n",
    "plt.figure()\n",
    "do_plot(1, RTM_TD, 'gray', 'Vp [m/s]', r\"Kleinneudorf p100 RTM result (TDFD DENISE: 40 - 50 Hz)\", -clip, clip)\n",
    "#do_plot(2, RTM_TD, 'gray', 'Vp [m/s]', r\" \", -clip, clip)\n",
    "plt.savefig('Kleinneudorf_RTM_DENISE.png', bbox_inches='tight', format='png', dpi=200)\n",
    "#plt.savefig('Kleinneudorf_RTM_DENISE.pdf', bbox_inches='tight', format='pdf')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
